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The dynamics of sharp interfaces separating two non-hydrostatically stressed solids 
is analyzed using the idea that the rate of mass transport across the interface is pro- 
portional to the thermodynamic potential difference across the interface. The solids 
are allowed to exchange mass by transforming one solid into the other, thermody- 
namic relations for the transformation of a mass element are derived and a linear 
stability analysis of the interface is carried out. The stability is shown to depend on 
the order of the phase transition occurring at the interface. Numerical simulations 
are performed in the non-linear regime to investigate the evolution and roughening of 
the interface. It is shown that even small contrasts in the referential densities of the 
solids may lead to the formation of finger like structures aligned with the principal 
direction of the far field stress. 

PACS numbers: 68.35.Ct, 68.35.Rh, 91.60.Hg 

I. INTRODUCTION 

The formation of complex patterns in stressed multiphase systems is a well known phe- 
nomenon. The important studies of Asaro and Tiller jl] and Grinfeld {^J brought attention 
to the morphological instability of stressed surfaces in contact with their melts or solutions. 
In the absence of surface tension, small perturbations of the surface increase in amplitude 
due to material diffusing along the surface from surface valleys, where the stress and chem- 
ical potential is high, to surrounding peaks where the stress and chemical potential is low. 
Important examples of instabilities at fluid-solid interfaces include defect nucleation and is- 



land growth in thin films 



4, so 



of fractal clusters by aggregation 



idification [5| and the formation of dendrites and growth 
The surface energy increases the chemical potential 
at regions of high curvature (convex with respect to the solution or melt, at the peaks) and 
reduces the chemical potential at region of low curvature (at the valleys) and this introduces 
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a characteristic scale below which the interface is stabilized. 

In systems where the fluid phase is replaced by another solid phase, i.e. solid-solid 
systems, the interface constraints alter the local equilibrium conditions. Here we study a 
general model for a propagating interface between non-hydrostatically stressed solids. The 
interface propagates by mass transformation from one phase into the other. The phase 
transformation is assumed to be local, i.e. the distance over which the solid is transported 
via surface diffusion or solvent mediated diffusion is negligible compared to other relevant 
scales of the system. Although the derivations apply to a diffuse interface, we shall here 
treat only coherent interfaces, where there is no nucleation of new phases or formation of 
gaps between the two solids 0,0], in the sharp interface limit. For example, in rocks such 
processes appear at the grain scale in "dry recrystallization" ji], lo|. Common examples of 
coherent interfaces that migrate under the influence of stress include the surfaces of coherent 
precipitates (stressed inclusion embedded in a crystal matrix) t| and interfaces associated 
with isochemical transformations. Most studies of solid-solid phase transformations have 
been limited to the calculation of chemical potentials in equilibrium and have provided 
little insight into the kinetics. Here we investigate the out of equilibrium dynamics of mass 
exchange between two distinct solid phases separated by a sharp interface. We expand on the 
recent work presented in 11] where we studied the phase transformation kinetics controlled 
by the Helmholtz free energy. It was shown that a morphological instability is triggered by 
a finite jump in the free energy density across the interface, and in the non-linear regime 
this leads to the formation of finger like structures aligned with the principal direction of 
the applied stress. 

In the majority of solid-solid phase transformation processes, the propagation of the 
interface is accompanied by a change in density. For this reason the density is an important 
order parameter that quantitatively characterizes the difference between the two phases. 
We consider two types of phase transitions underlying the kinetics, first order and second 
order, which result in fundamentally different behaviors at the phase boundary. A first 
order phase transition occurs when the two phases have different referential densities and 
it typically results in morphological instability along the boundary whereas a second order 
phase transition may either stabilize or destabilize the interface depending on Poisson's 
ratios of the two phases. A simple sketch of the stability diagram is outlined in Fig. [1] for 
relative values of density and shear modulus of the two phases. 
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FIG. 1: Sketch of a stability diagram for the growth rate of a sharp interface separating two solid 
materials. The axes show relative values of the shear modulus and density of the phases. As it will 
be shown in Sec. Ill, the symmetry of the diagram is broken by the values of the Poisson's ratios. 

The article consists of five sections. In Sec. II we derive a general equation for the 
kinetics for mass exchange at a solid-solid phase boundary separating two linear elastic 
solids. We utilize the derived equations on a simple one dimensional example and offer a short 
discussion of the order of the phase transition underlying the kinetics. We proceed in Sec. 
Ill with a linear stability analysis of the full two-dimensional problem. In two dimensions, 
the phase transformation kinetics gives rise to the development of complex patterns along 
the phase boundary. While we solve the problem analytically for small perturbations of a 
flat interface, things become more complicated in the non-linear regime, and we resort to 
numerical simulations based on the combination of a Galerkin finite element discretization 
with a level-set method for tracking the phase boundary. In Sec. IV, numerical results are 
presented together with discussions. Finally in Sec. V we offer concluding remarks. 

II. GENERAL PHASE TRANSFORMATION KINETICS 

Although the equations that we derive for the exchange of a mass element between two 
solid phases in a non-hydrostatically stressed system apply to more general settings, we limit 
ourselves to the study of two solids separated by a single sharp interface. The solids are 
stressed by an external uniaxial load as illustrated in Fig. |2j In the referential configuration, 
a solid phase is assumed to have a homogenous mass density, p°, defined per unit undeformed 
volume occupied by that phase. After the deformation, the densities are functions of space x 
and time t, i.e. pi(x, t) and p2^x, t). The average density of the two-phase system is denoted 
by p(x,t). Finally, the mass fraction for phase 1 is denoted by c. In this notation, the mass 
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fraction of phase 2 becomes 1 — c. 

For non-vanishing densities, the mass-averaged velocity is defined as 

v = cvx + (1 — c)v 2 . 



(1) 



Throughout the text, the mass average of any quantity is indicated by a bar. Similarly, the 
average specific free energy density is given by 



/=c/ 1 + (l-c)/ 2 . 



(2) 



The total specific volume is related to the real densities in the deformed state, p\(x,t) and 
p 2 (x,t) by 



P 



cpi 1 + (l-c)p 2 1 . 



(3) 



The interface separating the two phases is tracked by the zero level of a scalar field <f)(x, t) 
passively advected according to the equation 



dt 



+ W |V0| = 0, 



(4) 



where W is the normal velocity of the surface. It follows that the interface is given by the 
zero level set 

r = {x\<f>(x,t) = 0, for alU}. (5) 

The scalar field is constructed such that phase 1 occupies the domain in which <f)(x, t) > 
and phase 2 occupies the domain in which (p(x, t) < 0, see Fig. [2J In this notation, the mass 
fraction may be expressed as the characteristic function of the scalar field, 



c(x,t) =H((f)(x,t)) 



1, if 4>(x,t) > 
|, if <j>[x,t) = 
0, otherwise. 



(6) 



In the subsequent analysis, we make use of the following relations (see e.g. (l2j) 



VjC = Uidr, d t c = -W5 r , 



(7) 



where n« = Vj0/|V0| is the normal unit vector of the interface, W = — dt<f)/\ V0| is the 
normal velocity and 5-p = |V0|<5(0) is the surface delta function. 
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FIG. 2: (color online) Two solids separated by a sharp interface. A compressional force is applied 
at the margins in the vertical direction 

Taking the gradient of the averaged velocity from Eq. ([1]) and using the above identities, 
the following relation is obtained 



dv ■ 

ViVj = -Q^^i C + cVif ij + (1 - c)ViV 2t j 



dc 



riiSv + ViVj. 



(8) 



A. Kinetics of the phase transformation 

The system must satisfy fundamental conservation principles for the mass, momentum, 
energy and entropy. Let us denote the material time derivative with respect to the mass- 
averaged velocity by a dot, i.e. G = <9f0 + UjVj0. Then, the local mass conservation can be 
written in the form 

p = -pViVi. (9) 
and the local momentum balance can be written in the form 



pvi = Via. 



(10) 



where is the stress tensor. 
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The mass fraction of phase 1 satisfies the advection-reaction equation given by 

pc = Q5 r , (11) 

where the mass exchange rate Q is confined to the interface by the delta-function (in the 
sharp interface limit). Mass transport by diffusion is negligible in the reaction dominated 
regime. This is a valid approximation when the characteristic length I = T>/W, where T> is 
the diffusion coefficient and W is the velocity of the interface, is small compared with other 
relevant microscopic length scales. That is material diffusion occurs on a time scale much 
longer than any other relevant time scale in the system or equivalently the characteristic 
length scale formed from the diffusion constant and solidification or precipitation rate is 
small compared to other relevant microscopic scales. 

In the linear kinetics, the mass exchange rate is now derived from the requirement that the 
entropy production has a positive quadratic form. We start by expressing the conservation 
of specific energy density e in the form 

= (TijViVj, (12) 

where v 2 = cvf + (1 — c)v\ since the cross term vanishes in the limit of a sharp interface. 
At equilibrium 

e = f + Ts, (13) 

where the free energy is assumed to be a function of the local strain and the composition, 
i.e. / = f(eij,c). By inserting the energy conservation equation, Eq. (1121) . into the time 
derivative of this equation, under constant temperature conditions, the expression 

P TS = (TijViVj - Pq—^J ~ Pfc C i ( 14 ) 

is obtained. The phase transformation is assumed to be slow and isothermal. From Eqs. 
(12]) and (JSJ) it follows that 

Qy . Q f , Q f 

P T ^ = a "i-fc S r + °i3^i v 3 - Pg—^ ~ ~q c P L ( 15 ) 



Given that the strain rate is = l/2(ViVj + VjVi) and using the symmetry of the stress 
tensor, we arrive at the expression 
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where a n j = a^rii is the stress vector at the interface. From Eqs. (|HJ) and (Q and using an 

equation of state of the form pfej, c) = p°(c)(l — e&) it follows that, 

<9p . dp ^_ dv n — — 

Jf" c + J^ - = — fT p r ~ p iVi ^ 

o— — - 

-TrQ^r - p £ii = — ^— pc>r - pviVi => 
p oc oc 

Using Eq. (j3J) for the density, the jump in the material velocity is related to the reaction 
rate by 

" </:Tl (") 



9c 9c \p / 

The direction of the kinetics is constrained by the second law of thermodynamics which can 
be expressed in the continuum form as 

ps + vat = n„ (is) 

where J? is the entropy flux density and U s > is the entropy production rate. We consider 
the case where the entropy flux is negligible (in the absence of mass and heat fluxes) and 
therefore set J s = 0. Combining Eqs. f|T6|) and ffTSl . it can be seen that the positive entropy 
production rate leads to the condition 

*-£ G) - 1) Qh + (?« - p S w > = Tii - - ° (i9) 

on the reaction rate. We now define a constitutive relation that couples the stress to the 
strain via the Helmholtz free energy, 

. y =p|-. (20) 

From Eq. (1T9]) we observe that the entropy is produced only at the interface, and in the 
linear kinetics regime the reaction rate is proportional to (see e.g. 131]). 



dc \p J dc 
where K > is a system specific constant. 

The normal velocity of a sharp interface is obtained by integrating Eq. (jTTj) across the 
interface and taking the singular part of it, 



K 

W v n 

P 



l-f 

@nn J 
P 



(22) 
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Here we introduce the jump in the quantity a from one phase to another [a] := a\ — a 2 , 
where is the value of tij in phase i outside the interface zone as the interface is approached. 
The additional interfacial jump conditions of the total mass and force balance from Eqs. ([9]) 
and (ITUl) are given by 

lp(W-v n )}=0 (23) 
l* ijnj } = 0. (24) 

In general, surface energy 7 and surface stresses may have an important effect on the 
kinetics at the phase boundary with high curvature )C, therefore the expressions given above 
are modified to take this into account. For this purpose we utilize the Cahn-Hilliard for- 



malism 



14| of a diffuse interface. The surface energy is obtained by allowing the Helmholtz 



free energy density to be a function of the mass fraction gradients, i.e. 

Pf{<,r c, Vc) = pf (e, c) + y I Vc| 2 , (25) 

where K\ is a small parameter related to the infinitesimal thickness of the interface and fo 
is the homogenous free energy density introduced above. Because the composition gradient 
is small everywhere except for a thin zone at the interface, the free energy can be separated 
into bulk and surface contributions. If we now take the limit of vanishing surface thickness 
and follow the derivations in the appendix we obtain the general jump condition for the 
normal force vector, 

KnJ = -2/C7. (26) 

In the aforementioned expression of the interfacial velocity Eq. 0221) the normal stress vector 
was continuous across the interface. In the presence of surface tension, the normal velocity 
is altered by an additional contribution from the surface energy, 

W » v hn + - ([/] - W [p" 1 ] + 2/C 7 (p- 1 » , (27) 
Pi 

where we have used the interface average defined as (a) = l/2(ai + a 2 ). 

B. Example: Phase transformation kinetics in a one dimensional system 



We start out considering the phase transformation kinetics of a one dimensional system 
composed of two linear elastic solids separated by a single interface. A force a is applied 



E 2 , L 2 , m 2 



G 



FIG. 3: One dimensional system undergoing phase transformation 



at the boundary of the system (see Fig. [3]) and each solid phase is represented by its 
Young's modulus Ei {i = 1, 2), undeformed density p° and length Vj. In the deformed state 
when the external force is applied the length becomes Lj = L° (1 + a/E^) and the density 
Pi = p^Vf/Li. The specific free energy is given by 

.2 / 1 _ £ . 



f = °-( c - + 

2 \p 1 (E 1 +<r) p 2 (E 2 + a) 



(28) 



In the following, we do not allow new phases to nucleate within the solids and limit our 
considerations to the propagation of a single interface separating the solids. The system is 
assumed to be isothermal and no diffusion of mass takes place. The interface moves as one 
phase, slowly transforms into the other and an amount p\dLi, of solid 1 is replaced by an 
amount pidL^ of solid 2, with conservation of the total mass. The phase transformation is 
assumed to be irreversible and to occur on time scales that are much larger than the time 
it takes for the system to relax mechanically under the deformational stresses. 

In the one dimensional setting the local mass exchange rate is given by a linear kinetic 
equation, Eq. (12~T1) . of the form 



mi 



-K 



[ a 2 a 




'a 2 a j 


= K 


[2p°E p 


2p°E + p°J 



(29) 



with K > 0. In most cases, the contribution from the jump in the elastic energy density 
will be small compared to the contribution from the work term (because a/E 1, within 
the linear elasticity regime). The change in the total length will in general follow the sign 
of the stress 



Li(l 
K 



Pi- 
Pi 



mi 



a- 



a 



2Ep° + p 



a 



IP 



+ Ep° 



If the densities in the undeformed states are identical, p\ = p 2 , the change in the total length 
is given by 



K 



2> 



1 

E 



(30) 
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whereas a jump in the referential densities [p\ 7^ p°) wm result in a work term given by 

(31) 

Under a compressional load, the dense phase grows at the expense of the less dense phase 
(if the two phases have the same Young's modulus) and the soft phase grows at the expense 
of the hard phase (if the two phases have the same density), such that overall the system 
responds to the external force by shrinking. The one-dimensional model cannot predict the 
morphological stability of the propagating phase boundary in two dimensions. It turns out 
that the work term destabilizes the propagating boundary under a compressional load. 



C. First and second order phase transitions: Equilibrium phase diagrams 

In the above derivations, the reaction rate is determined by the jump in the Gibbs poten- 
tial across the phase boundary. Whenever the system is stressed, only one of the two phases 
will be stable, i.e. the general two phase system will always evolve to an equilibrium state 
consisting of a single phase. In the absence of an external stress, it is possible for two phases 
to coexist without any phase transformation taking place. In the one dimensional example, 
the relevant field variable is the stress o applied to the system and the Gibbs potential is 
given by (follows from Eq. (129]) ) 

g(a) = - -. (32) 

In Fig. H]we show an equilibrium phase diagram in the conjugate pair of variables a and 1/p. 
If the derivative of the Gibbs potential with respect to the external field a is evaluated at the 
critical point a = 0, it can be seen that there are two possible scenarios. The first scenario 
is a first order phase transition, which occurs whenever there is a jump in the referential 
densities, i.e. the derivative of the Gibbs potential is discontinuous and the second derivative 
diverges at the critical point. The other scenario is a second order phase transition, which 
occurs when the referential densities of the two phases are identical. We then have a jump 
in the second order derivative whenever Young's modules of the two phases are dissimilar. 

The order of the phase transition has a fundamental impact on the dynamics. In two 
dimensions a first order phase transition kinetics will generally lead to morphological insta- 
bilities of the propagating phase boundary while a second order phase transition will either 
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FIG. 4: (color online) Part (a) illustrates the phase diagram for a second order phase transition in 
the p — a plane. The solid-solid kinetics will always be directed from the unstable phase (dashed 
line) to the stable phase as illustrated by reaction path Q marked by the dashed arrow. The slopes 
of the densities with respect to stress are Young's modules of the materials. Part (b) illustrates the 
equilibrium curves of the first order phase transition. For the first order phase transition one would 
in general expect to see hysteresis effects extending the continuous lines (stable regions) beyond 
the point a = 0. Here we have shown an idealized case where such effects are disregarded. 

flatten or roughen the boundary depending on Poisson's ratios of the two materials. In 
the next section we analyze the different phase transitions by performing a linear stability 
analysis. 

III. LINEAR PERTURBATION ANALYSIS 

We now solve the elasto-static Eqs. ffTUl) and fliZoT) together with the kinetics Eqs. f[2"2"j) 
and (12T1) in two dimensions for an arbitrary perturbation to an initially flat interface using 
the quasi-static version of momentum balance in Eq. (jTUI) . In addition to the translational 
dynamics observed in the one-dimensional system presented above, it turns out, that in two 
dimensions the interface dynamics is non-trivial and may lead to the formation of finger-like 
structures. The general setup is shown in Fig. [2] where phase i, i — 1, 2, has material 
parameters /Xj, z/j and pi, with /ij being the shear modulus and z/j being the Poisson's ratio. 
In general, the interface velocity depends on its morphology, the 6 material parameters and 
the external loading a^. One degree of freedom is removed by rescaling the shear modulus 
of one phase with the external load. 
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A. Stress field around a perturbed flat interface 

In order to evaluate the jump in Gibbs energy density, i.e. + W], we need to 

determine the stress field around the interface by solving the elastostatic equations. We 
have that under plane stress conditions, the local strain energy density can be written on 
the form 

T = 4^ + a yy ~ \T^ {axx + ayy)2 + 2a *y) (33) 

and the work term is defined as 

W = -(Jnnpr 1 = -<r nn p$(l + Tr(e)). (34) 
The trace of strain is given in terms of stress by 

Tr ( e ) = 2^(1+1) ( Cxx + ° yy) - ^ 

Note that we could as well have formulated the problem under plane strain conditions; 
however, the generic behavior in both plane stress and strain is the same although the 
detailed dependence on the material parameters is altered. 

We solve the mechanical problem by finding the Airy stress function, U(x,y) [15], which 
satisfies the biharmonic equation A 2 U = 0. Once the stress function has been found, the 
stress tensor components readily follow from the relations 

&u_ mj_ du 

<Txx- dy2 , °yy ^ , °xy ^ W 

The biharmonic equation is solved under the boundary conditions of a normal load applied 
in the y direction at infinity, i.e. a yy — > —\<Joo\ < and a xy = for y — > ±oo. The continuity 
of the stress vector across the interface follows from force balance. In addition we require 
that u x (±oa,y) = 0. 

For a flat interface, the stress field is homogenous in space. This implies that the Airy 
stress function is quadratic in x and y, with coefficients determined by the boundary condi- 
tions. With the boundary conditions specified above, the stress function for the i-th phase 
can be written in the form 

U l (x,y)= l -^(x 2 + v l y 2 ). (37) 

From this stress function we can calculate the Gibbs potential which in the case of dissim- 
ilar phases is discontinuous across the interface. The velocity of the phase transformation 
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readily follows from the potential 



W oc w; + wo]=ki(44) 

VPl H2/ 



4 V PiPi P2P2 

The subscript of the free energy density and the work term refers to an unperturbed interface. 
From the above equation, we see that when the lower phase is much denser than the upper 
phase, i.e. p\ <C p\, the interface propagates uniformly into the upper phase with a velocity 
W ~ | (Too | [1/p] > 0, i.e. the denser phase grows into the softer. When the densities are 
identical or almost identical, p 2 /pi ~ 1 and the shear modules significantly different, i.e. 
Pi P2- When the two solids phases have identical Poisson's ratios, z/, we see that the 
softer phase can only grow into the harder one when v < 1/3. 

In the case of an arbitrarily shaped interface separating the two phases, the analytical 
solution to the stress field is in general far from trivial. In-plane problems can in some cases 
be solves using conformal mappings or perturbation schemes Ha, ll6|, Il7| . Here, we solve 
the stress field around a small undulation of flat interface employing a linear perturbation 
scheme 17J. In the linear stability analysis we now study the growth of an arbitrary harmonic 
perturbation with wavelength k, i.e. h(x,t) = Ae ut cos(kx) with A <C 1. In appendix [Bj we 
derive expressions for a general perturbation. The Airy stress function can be written as a 
superposition of the solution to the flat interface and a small correction due to undulation, 
U(x, y) = Uq(x, y) + Q(x, y), where 0(x, y) is determined from the interfacial constraints of 
continuous stress vector and displacement field. When the wave number k is much smaller 
than the cutoff introduced by the surface tension, we obtain the following expressions for 
the Airy stress functions 

\aoo\hix) exp(-%)(aiy + (3) 



k{\x 2 K\ + pi)(pi«2 + P2) 



(Too \h(x) exp(ky)(a 2 y - fl) 
®2{x,y) = — ■ ■ — (39) 

fc(//2«l + Pl)(Pl«2 + P2j 

where Ki = and we have introduced the material specific constants, 



a x = k(l - za)(P2 - ni)(piK2 + p 2 ) 
a 2 = k(l - V2)(nx - £t2)(/42«i 
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and 



2f4 



! ' 2 ~2f4 



l + Vx 



v x - v-i 



1 + 1/2) (1 + ^1) 



From the Airy stress functions, we then calculate the stress components using Eq. fl36|) 
and find the jumps in the Gibbs energy density from Eqs. (f53|) and (1541) . The evolution of 
the shape perturbation relative to a uniform translation of the flat interface is described by 
Eq. ff271) . namely 

^|^oc[^ + W]-^o, (40) 
which in the linear regime corresponds to a dispersion relation given in the general form as 

Ul OC ; . (41) 

h 

Below follows an evaluation of the growth rate for a small harmonic perturbation to a flat 
interface. For this perturbation, the general expression for the growth rate follows directly 
upon insertion of the Airy functions in Eq. ( 1391) and then in Eq. ( |36i) . however, the growth 
rate is not easily expressed in a short and readable form and we have therefore limited 
our presentation to a few special cases. The growth rate is a function of the six material 
parameters [y^ p i: pi) and the external stress. Naturally, the stability of the growing interface 
is invariant under the interchange of the solid phases and correspondingly the region of the 
stability diagram that we have to study is reduced. 



B. First and second order phase transition: Stability diagrams 

In the second order phase transition when both solids have the same referential densities 
Pi = P2 = P° and when the Poisson's ratios v\ — v-i — v are identical the dispersion relation 
assumes a simple form 

u = (3z/ - 1)(1 - v){pi + p 2 ){p2 - Pi) 2 ^ 
k pViA*2(A*i+/i2K)(A*2 + A*i«)(l + ^) 

where k is the fraction introduced above and k the wave number of the perturbation. The 

expression reveals an interesting behavior where the interface is stable for Poisson's ratio 

less than 1/3 and is unstable for Poisson's ratio larger than 1/3. Fig. (JSj) shows stability 

diagrams for the specific case where pi = 1 and p\ = 1 (in arbitrary units). In panel (A) 
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the diagram is calculated for two solids that have the same Poisson's ratio and with a value 
v = 1/4. The second order phase transition occurs along the horizontal cut p° = 1 and 
is marked by a dashed grey line. We observe that u/k is negative along this line and the 
interface is therefore stable. For v larger than 1/3 (not shown in the figure) the horizontal 
zero level curve will flip around and the grey dashed line will then be covered with unstable 
regions. In order to see this flip, we expand Eq. (j4Tj) around the point (1,1), i.e. in terms of 
p\ — 1 and /i 2 — 1, and achieve the following expression for the zero curve 



„ , , (l-2^-3^ 2 )( M2 -l) 



Pi ~ 1 + 7^ -T^ (43) 

2 u(7 + u) 

Note that the right hand side is in units of p\. We directly observe that the horizontal 
zero curve flips around at the critical point v = 1/3. In the case when the two solids are 
identical, i.e. at the point (1,1) in the stability diagram, all modes will as expected remain 
unchanged and the interface therefore remain unaltered. The other parts of the zero levels 
lead to marginal stability but will in general induce a growth of the interface due to the 
unperturbed Gibbs potential Eq. fl38l) . We now consider a cut in the stability diagram 
where the two solids have the same shear modules, p.\ — p-i — A*, but different densities and 
Poisson's ratios. For different Poisson's ratios the dispersion relation Eq. (|4ip becomes 

^ = (^2 - v\){v\p\ ~ V2P1 + 2(p 2 l - p\)p) 

k Ap\p»p 1 > 

From this expression we see that the vertical zero line observed in Eq. (j42H and in Fig. [5] 
panel (A) only exists for identical Poisson's ratios. When the solids have different Poisson's 
ratios, the separatrix or intersection of the two zero curves located at (1,1) in panel (A) will 
split into two non-intersecting zero curves. In panel (B) we show a stability diagram for 
solids with Poisson's ratios V\ = 0.45 and V2 = 0.40. 

In general the stability diagram is characterized by four quadrants, two stable and two 
unstable, delimited by neutral zero curves. The physical regions would typically correspond 
to the quadrants / and III under the assumption that higher density implies higher shear 
modulus. In these quadrants the growth rate is typically positive (i.e. the interface is 
unstable) except for a thin region at the borderline between a first and second order phase 
transition, i.e. when p 2 — p%. 



16 



Q. 1 




1 1.5 
H2, (^1 = 1) 



0.05 



0.00 



-0.15 



Q_ 1 - 



0.5 



J I I C\l c 

jL(§U cj> <J gf c 

r> ?' ?• 


O O (2, 
N> > <3> 






°o<> 














0.5 



1 1.5 

M-2. (M-1 = 1 ) 



0.20 



0.15 



-0.05 



-0.10 



FIG. 5: (color online) Panel (A), stability diagram for two solids materials with identical Poisson's 
ratio of v = 0.25. Panel (B), diagram for solids with Poisson's ratios of v\ = 0.45 and = 0.40. 

IV. NUMERICAL RESULTS AND DISCUSSIONS 



The linear stability analysis revealed an intricate change in stability depending on the 
material properties and densities of the two solids. We explore this stability beyond the 
linear regime using numerical methods. The bulk elasto-static equation Eq. tfTOl) is solved 
numerically on an unstructured triangular grid using the Galerkin finite element method 
and the surface tension force is converted to a body force in a narrow band surrounding 
the interface. The discontinuous jumps appearing in the dynamical Eq. (|27j) are computed 
at the outer border of the band. Periodic boundary conditions are used to minimize the 
possible influence of the finite system size in the x-direction (parallel to the interface). The 
interface is tracked using a level set method (e.g. isj]) and propagated with the normal 
velocity calculated in Sec. II using Eq. (|27j) . Several level set functions, <p(x, t), can be used, 
however, most level set methods use the signed distance function (\<f>(x, t)\ is the shortest 
distance between x and the interface and the sign of <fi(x, t) identifies the phase at position x). 
Good numerical accuracy can be obtained by keeping 4>{x, t) a signed distance function at 
all times, and this is achieved by frequent reinitialization of 4>(x, t) according to the iterative 
scheme 

^ + 5(0 O )(|V0|-1) = O, (45) 
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where 0o is the level set function before the reinitialization, t' is a fictitious time, and 



are needed at each time step, to obtain a good approximation to a signed distance function, 
and it is only necessary to update the level set function in a narrow band around the interface. 

In Figs. ([6]) and (J7J) we present numerical simulations of the phase transformation kinetics 
using parameter regions where the interface is either stable or unstable. The simulations 
presented in Fig. ([6]) panels (A) and (B) represent interface snap shots of a first order phase 
transition dynamics and panels (C) and (D) simulations of a second order, respectively. In 
panel (A), the values of the parameters were chosen in a region of the stability diagram 
where the interface is predicted to roughen and in panel (B) we have used parameters 
corresponding to a stable evolution of the interface. Note that the interface in both cases 
is moving from the dense phase into the soft phase independent of its stability. This is in 
agreement with the one dimensional calculation performed in Sec. II. Panel (C) shows a case 
of a second order phase transition where the interface is unstable, while panel (D) shows 
a stable case. We notice that, for second order phase transitions, the overall translation of 
the interface is changed in unison with its stability. In Eq. ( )43l) we saw that the stability of 
the second order phase transition is dictated by the values of Poisson's ratios. For Poisson's 
ratio smaller than 1/3, the kinetics is stable and the phase of small shear modulus grows into 
the phase of higher shear modulus while for higher values of Poisson's ratio the behavior is 
reversed and the interface roughens with time. This also follows from Eq. (1381) . In fig. we 
have plotted the mean velocity as a function of time for the simulations presented in fig. El 



In conclusion, it has been shown that the phase transformation of one solid into the 
another across a thin interface may lead to a morphological instability, as well as the devel- 
opment of fingers along the propagating interface. We have presented a stability analysis 
based on the Gibbs potential for non-hydrostatically stressed solids and have established 
a linear relationship between the rate of entropy production at the interface and the rate 
of mass exchange between the solid phases. The solids are compressed transverse to the 
interface and corresponding stability diagrams reveal an intricate dependence of the sta- 
bility on the material density, Poisson's ratio and Young's modulus. With the density as 




Generally only a couple of iterations 



V. CONCLUDING REMARKS 
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FIG. 6: (color online) Simulations of the temporal evolution of solid-solid interfaces for first order 
(Panels (A) and (B)) and second order (Panels (C) and (D)) phase transitions. Panel (A) shows a 
simulation using p\ = 1.0, /xi = 1.0 and p2 = 1.05, p,2 = 2.0. Both phases have identical Poisson's 
ratio, ui = V2 = 0.45. Panel (B) is a simulation run with densities and shear modules similar 
to panel (A) but with a different Poisson's ratios, v\ = vi = 0.25. Panel (C) is a simulation 
run with p\ = 1.0, /xi = 1.0 and p2 = 1-0,^2 = 2.0. Both phases have identical Poisson's ratios, 
v\ = V2 = 0.45. Panel (D) shows a simulation run with densities and shear modules similar to 
Panel (C) but with different Poisson's ratios, V\ = v% = 0.25. The color code represents a time 
arrow pointing from the darker blue (early stage) to the lighter blue (final stage). 

order parameter, two types of phase transitions were considered, a first and a second order, 
respectively. 

For both types of transitions we find expressions for the curves separating the stable 
and unstable regions in the stability diagram. For most material parameters the first order 
phase transition, i.e. when the two solids have different referential densities, destabilizes 
the interface by allowing fingers to grow from the denser phase into the other. When the 
solids have identical or almost identical densities, i.e. a second order phase transition, we 
find that the stability depends on Poisson's ratios of the two solids. If the two solids have 
Poisson's ratios less than 1/3, the phase transition dynamics of the two solids will lead to a 
flattening of the interface, i.e. any perturbation of a flat interface will decay and ultimately 
the interface will propagate uniformly from the soft phase (low Young's modulus) into the 
hard phase (high Young's modulus). 
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FIG. 7: (color online) Normal velocity as a function of time in first order (Panels (A) and (B)) and 
second order (Panels (C) and (D)) phase transitions. The simulations presented in the individual 
panels are identical to the corresponding panels in Fig. The color code of the curves reference 
to the mean velocity (in black), the mean of lowest 10 % (in red) and the mean of the highest 10 
% (in green). 

We believe that our classification of the phase transition order together with the stability 
analysis may find application in many natural systems, since the morphological stability 
directly provide information about the order of the underlying phase transformation process 
and the material parameters. 
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APPENDIX A: SURFACE TENSION 

In this appendix we present additional details on the derivation of the reaction rate Eq. 
(I27p including the interfacial free energy. Let us consider a diffuse interface characterized by 
a small thickness over which the concentration field varies smoothly between the constant 
values in the bulk of the two phases. In the Cahn-Hilliard formalism, the free energy is 
introduced as a function of both the concentration and concentration gradients, and has the 
form 

p/(ey,c, Vc) = p/o(e y ,c) + y|Vc| 2 , (Al) 

where the first term is the free energy in the bulk and the second term is associated with the 
interfacial free energy. Here K\ is a small parameter related to the thickness of the interface. 

In this case, the calculation of the reaction rate Q proceeds as in Sec. II. We apply the 
total time derivative of the local equilibrium equation, Eq. (Tl3l) . where the free energy is 
given by Eq. (IA1I) and then obtain the following expression 

df df df 
e = + -fc + — f- V;c - VjcVtVj) + Ts, (A2) 

den oc oViC 



where the commutation relation, 4VjC = VjC — V{V jV^c, has been used 



121 ] . Combining 



the above equation with the conservation of energy from Eq. (1121) and the entropy balance 
from Eq.( fl8|) an expression for the entropy production rate is obtained 



, df \ fdf df \ . df _ 



= m [a l3 + P V t c— ) n 3 QS r - ) Q5r 

( „ df df\^_ 

We observe that Il s satisfies the second law of thermodynamics provided that the last term 
vanishes and the rest of the terms are brought into a quadratic form. This implies a consti- 
tutive equation for the stress given by 

9f df / A q\ 

a - =P de- 3 - pVlC dV7c (A3) 
and a linear kinetics law with the reaction rate being proportional to 

n~v( df° 9 fl\ df df \ 
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where K is a positive local constant of proportionality and is the elastic stress in the 
absence of surface tension. 

Using Eq. (ED), the two constitutive laws may be expressed as 

Oij = (T°- - KjViC ® VjC (A5) 

« - * 6) - f + • < A6 > 

where of,- is the elastic stress obtained in Sec. II without the surface stress. 

In the sharp interface limit, i.e. the thickness goes to zero, the surface free energy becomes 

p/ -/ = Kl |Vc| 2 ^75r, (A7) 

and surface stress is related to the surface energy by 

= -ilVcl 2 (l - ^ ^) - 7(1 - m® n 3 )5 T . (A8) 

The divergence of the surface stress is then calculated as 

V,^ r/ = 2/C 7 Vr, (A9) 

where /C is the local curvature. 



APPENDIX B: GOURSAT FUNCTIONS AROUND A PERTURBED FLAT 

INTERFACE 

In this appendix, we explain in details how to calculate the Airy stress functions around 
the perturbed flat interface introduced in Sec. III. All the detailed calculations were carried 
out in Maple in order to handle the lengthy algebraic expressions. 

The Airy stress function satisfies the biharmonic equation d\d\V = 0. This equation has a 
general solution which can be written in the Goursat form U(z, z) = 3l{z<f)(z) +x( z )}: where 
<p(z) and x{ z ) are complex functions determined by the boundary conditions. Combining 
Eq. (|36|) with the Goursat solution, stress components are related to these functions by the 
following expressions 

cr(z) = a xx (x,y) + ayy(x,y) = 2{ip'(z) +ip'(z)}, (Bl) 
E(z) = cryy(x, y) - a xx (x, y) + 2ia xy (x, y) 

= 2{zcp"{z)+TP{z)}, (B2) 
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where f{z) = x'{ z )- The solution to the biharmonic equation is determined up to a linear 
gauge transformation, 

<p(z) (-»• <p(z) + Ciz + p (B3) 
ip(z) ^ ip(z) + q, (B4) 

where C is a real number and p, q are arbitrary complex numbers. 

The boundary conditions are given by the far-field stresses and the constraints at the 
interface. Here we consider that the system is loaded by a uniaxial compression in the y- 
direction, a yy (x,oo) = — [a^l < 0. Whenever the two phases are different an interface is 
introduced at which we require force balance and continuous displacement field. The force 
balance is expressed by the following jump condition 

i<J xx n x + a xy n y + i(a yx n x + cr yy n y )} = -jfC(n x + in y ), 

where JC is the local curvature and 7 is the surface tension. From Eqs. ( 1B1I) and ( 1B2I) we 
find that the force balance leads to the following condition on the Goursat functions 

[</? + zip' + vpj=i 7/C(n x + in y )ds, (B5) 
Jo 

where s is a point at the interface. The continuity of the displacement field across the 
interface introduces an additional jump condition given by 



-(-Kip + Zip' + *0) 



0, (B6) 



where \x is the shear modulus and constant for in-plane stress-elasticity deter- 

mined by the Poisson's ratio. 

The two jump conditions, Eqs. (JB5]) and (1B6P combined with the far-field boundary 
conditions, ipca(z) = — j('i- + i / )\o'oo\z and ipoo(z) = — — ^Icool-z are sufficient to determine 
the fields ipi(z), ipi(z), ¥2(2) and ^2(2). 

Superimposing an arbitrary perturbation with amplitude h(x) on the flat interface, the 
Goursat functions are slightly altered. They can be expanded to linear order in h(x) as 



follows [171 ]. 



ip(x) ta ip (x) + ih(x)tp' (x) + (B7) 
ip(x) « ipo(x) + ih(x)ip' (x) + (B8) 

(B9) 
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and *f>(x) are functions of h(x). Inserting this expansion into Eqs. (1B6j) and (1B5[) . we 
obtain that the corresponding jump conditions for the perturbation fields 



(BIO) 
(Bll) 



$0) + x$'(x) + *(x)J = ih(x) [S (x 

+ /(*) 
= ih(x 



-k$(x) + x$'(a;) + *ff(x) 



where f(x) = i jJC(n x + in y )ds. To linear order we find that f(x) ~ —7 f£ h"(s)ds. Eqs. 
(IBlOjl and (1B11I) can be rewritten equivalently as 



$i(x) -O + - (1 + A)$ 2 (x 



-ifi/i(x)S i(x) + * + / (x) 

1 + K 



$ 2 (x) - n(^x$ 2 (x) + # 2 (x)^ - (1 + A)$!(x) 



-«n/i(x)s 02 (x) - 1+ f(x). 

1 + K 



(B12) 



(B13) 



The constants appearing above are expressed in terms of the elastic moduli. Adopting the 
notation of [l^], these are given by 



A = K 

A = K 



1/H2 - 

\j\x 2 + k/hi 

- \j\x 2 



n 



n 



i//i 2 - 1/jJd 

k)\x 2 + l//ii 

- 1//J.2 



(B14) 
(B15) 



Eqs. (IB12|) and flB13j) are solved at an arbitrary point 21 in the complex plane by applying 
the Cauchy integral and using the analytic continuation of each function [15j . Let us denote 
the Cauchy integral over the perturbation amplitude 



Hi(z) 
H2(z) 



J_ f 3?Ldx, with 3?(z) > 
2ni J x — z 



h(x) 



■dx, with Q(z) < 0. 



2iri J x — z 

Notice that the two functions satisfy the following relations 



(B16) 
(B17) 



H 1 (z) = -H 2 (z), H 2 {z) = -H l {z) 

where the principal value of the Cauchy integral is considered when x is a point on the real 
axis. 
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Thus, by applying the Cauchy integral with $s(z) > in Eq. IB12I and $t(z) < in Eq. 
IB13t $i and ty 2 are determined in the integral form as follows 

= -in^ QA H 1 (z) + \^F 1 (z) 

1 ~~\- K 

$ 2 (z) = ^nS , 2 F 2 (z) + \±^F 2 (z), 



where 

1 ^ 2 



F'(z) = / ^^dx « -7— #(A (B18) 

^1(2) is calculated from the complex conjugation of Eq. ( 1B12I) when the Cauchy integral is 
applied on both sides of the equation and 3 (z) > 0. In a similar manner, $ 2 i z ) is derived 
from Eq. flB13j) . The final expressions for the two functions then follow 



#1(3) = -i^ Q , 1 H 1 (z)-\^(-iUE , 2 H 1 (z) 



1 + K ' 7 0(1 + k) 

* 2 {z) = iL 0fl H 2 {z)-^-^(inE 0tl H 2 {z) 



i + K v 7 n(i + «) 

For a cosine perturbation of the interface, h(x) = Acos(kx), with A <C 1 the Airy stress 
function, U(x,y) = ?R,{z(j)(z) + x( z )} is obtained explicitly. 
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